One-way analysis of variance: Multiple comparisons

DATAX121-23A (HAM) & (SEC) - Introduction to Statistical Methods

Briefly: The F-test’s for CS 2.2

Recall that

\[ \begin{aligned} H_0\!: & ~ \mu_\text{Low} = \mu_\text{Medium} = \mu_\text{High} \hspace{2.25em} \\ H_1\!: & ~ \text{At least one} ~ \mu_i \neq \mu_j \end{aligned} \]

We had strong evidence against the null that … in favour of the alternative that … (p-value = 0.01208)

# Fit the means-only model to the data
lm(GillRate ~ Calcium, data = respiration.df) |>
  # Decompose the total "variability" between and within groups
  anova()

Analysis of Variance Table

Response: GillRate
          Df  Sum Sq Mean Sq F value  Pr(>F)  
Calcium    2  2037.2 1018.61  4.6484 0.01208 *
Residuals 87 19064.3  219.13                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F-test ⇝ Multiple comparisons

Next is to construct confidence intervals for each pair of \(\mu_i - \mu_j\)
However, we cannot default to the same method taught in T04B¹. Why?
The goal is to construct a set of (1 - α)% confidence intervals that “match” the result of the F-test

Adjusting for multiple comparisons

To use ALL the data to quantify the uncertainty of each group’s sample mean, \(\text{se}(\bar{x}_i)\), we need to meet this assumption in particular

All groups have a similar measure of spread

Conveniently, the statistic we require was already calculated for the F-test: The mean square for residuals, \(MSR\)

\[ MSR = \frac{SSR}{n-k} \]

The intuitive rationale for MSR

Consider a randomised experiment with mice where a numeric response variable is of interest and say we split them into four treatment groups

Definition: (1 - α)% Confidence interval for μ_i − μ_j

\[ \bar{x}_i - \bar{x}_j \pm \frac{\text{Tukey}^*_{1-\alpha}(k, \nu)}{\sqrt{2}} \times \text{se}_{MSR}(\bar{x}_i - \bar{x}_j) \]

where:

\(\bar{x}_i ~ \& ~ \bar{x}_j\) are the sample means of the ith and jth group, respectively
\(n_i ~ \& ~ n_j\) are the number of observations in the ith and jth group, respectively
The confidence level is \((1 - \alpha)\), where \(\alpha\) is a proportion
\(\text{Tukey}^*_{1-\alpha}(k, \nu)\) is the Tukey multiplier for the prescribed confidence level of \((1 - \alpha)\)
- The number of groups, \(k\), that is, the number of levels for the categorical explanatory variable
- The degrees of freedom, \(\nu = n - k\)
\(\text{se}_{MSR}(\bar{x}_i - \bar{x}_j)\) is the standard error of \(\bar{x}_i - \bar{x}_j\) based on the mean square for residuals
- \(\text{se}_{MSR}(\bar{x}_i - \bar{x}_j) = \sqrt{MSR \times \left(\frac{1}{n_i} + \frac{1}{n_j}\right)}\)

CS 2.2 revisited: Fish respiration rates

# For the second time in DATAX121, we will load another R package
library(emmeans)

# Fit the means-only model to the data
respiration.fit <- lm(GillRate ~ Calcium, data = respiration.df)

# The emmeans package requires us to state categorical explanatory variable
emmeans(respiration.fit, ~ Calcium) |>
  pairs(infer = TRUE)

 contrast      estimate   SE df lower.CL upper.CL t.ratio p.value
 High - Low      -10.33 3.82 87   -19.45    -1.22  -2.704  0.0223
 High - Medium    -0.50 3.82 87    -9.61     8.61  -0.131  0.9906
 Low - Medium      9.83 3.82 87     0.72    18.95   2.573  0.0313

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 3 estimates 
P value adjustment: tukey method for comparing a family of 3 estimates

If this line causes an error, see Workshop 9 (when available)

The emmeans package expects a lm() object, which is why we now “saved” it output to a R object named respiration.fit

The contrast column of the R output tells us the order of the pairwise comparison

Interpretation of a confidence interval for μ_i − μ_j

# The emmeans package requires us to state categorical explanatory variable
emmeans(respiration.fit, ~ Calcium) |>
  pairs(infer = TRUE)

 contrast      estimate   SE df lower.CL upper.CL t.ratio p.value
 High - Low      -10.33 3.82 87   -19.45    -1.22  -2.704  0.0223
 High - Medium    -0.50 3.82 87    -9.61     8.61  -0.131  0.9906
 Low - Medium      9.83 3.82 87     0.72    18.95   2.573  0.0313

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 3 estimates 
P value adjustment: tukey method for comparing a family of 3 estimates

Rank order of the observed group means

Another question we could answer with the confidence intervals is to infer if there is any group(s) that were significantly different from the rest

Recall the output was

# The emmeans package requires us to state categorical explanatory variable
emmeans(respiration.fit, ~ Calcium) |>
  pairs(infer = TRUE)

 contrast      estimate   SE df lower.CL upper.CL t.ratio p.value
 High - Low      -10.33 3.82 87   -19.45    -1.22  -2.704  0.0223
 High - Medium    -0.50 3.82 87    -9.61     8.61  -0.131  0.9906
 Low - Medium      9.83 3.82 87     0.72    18.95   2.573  0.0313

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 3 estimates 
P value adjustment: tukey method for comparing a family of 3 estimates

Briefly: The adjusted p-value

The hypothesis test for each pairwise comparison is: \[ \begin{aligned} H_0\!: \mu_i - \mu_j = 0 \\ H_1\!: \mu_i - \mu_j \neq 0 \end{aligned} \]
The calculation for the set of p-values also need to “match” the result of the F-test¹
Like the confidence intervals, the calculation of the test statistic and the probability statement that calculates the p-value makes use of ALL the data

Interpretation of an adjusted p-value

# The emmeans package requires us to state categorical explanatory variable
emmeans(respiration.fit, ~ Calcium) |>
  pairs(infer = TRUE)

 contrast      estimate   SE df lower.CL upper.CL t.ratio p.value
 High - Low      -10.33 3.82 87   -19.45    -1.22  -2.704  0.0223
 High - Medium    -0.50 3.82 87    -9.61     8.61  -0.131  0.9906
 Low - Medium      9.83 3.82 87     0.72    18.95   2.573  0.0313

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 3 estimates 
P value adjustment: tukey method for comparing a family of 3 estimates